function W = lambertw3D(obj, beamshape)
% beamshape --- normalized in the incident xy-plane, dimensionless.
% res --- polarization distribution of the pumping field, dimensionless.

% normalized beamshape
funcVal = obj.getFunctionValue(beamshape);
nFac = simp3D(obj,funcVal);
funcVal = funcVal/nFac;

alpha0 = obj.physPara.power*funcVal/obj.physPara.P0;alpha0(alpha0>500)=0;
alpha1 = obj.physPara.power*funcVal/obj.physPara.P0;alpha1(alpha1<=500)=0;

W = (alpha1 + obj.physPara.OD*obj.unitGrid.Z.*alpha1./(1-alpha1))+...
     lambertw(alpha0.*exp(alpha0-obj.physPara.OD*obj.unitGrid.Z));
end

